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A  method  is  evaluated  for  estimating  the  absorption  coefficient  a  and  the  backscattering  coefficient  bb 
from  measurements  of  the  upward  and  downward  irradiances  Eu(z)  and  Ed(z).  With  this  method,  the 
reflectance  ratio  R(z)  and  the  downward  diffuse  attenuation  coefficient  Kd(z)  obtained  from  Eu(z )  and 
Ed(z )  are  used  to  estimate  the  inherent  optical  properties  and  Kx  that  are  the  asymptotic  values  of 
R(z)  and  Kd(z),  respectively.  For  an  assumed  scattering  phase  function  |3,  there  are  unique  correlations 
between  the  values  of  Rx  and  Km  and  those  of  a  and  bh  that  can  be  derived  from  the  radiative  transfer 
equation.  Good  estimates  of  a  and  the  Gordon  parameter  G  =  bb/(a  +  bb)  can  be  obtained  from  Rx  and 
K,r  if  the  true  scattering  phase  function  is  not  greatly  different  from  the  assumed  function.  The  method 
works  best  in  deep,  homogeneous  waters,  but  can  be  applied  to  some  cases  of  stratified  waters.  To 
improve  performance  in  shallow  waters  where  bottom  effects  are  important,  the  deep-  and  shallow- 
measurement  reflectance  models  also  are  developed.  ©  1997  Optical  Society  of  America 
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1.  Introduction 

Determination  of  the  beam  absorption,  scattering, 
and  backscattering  coefficients  a,  b,  and  bh  of  natural 
waters  is  a  primary  goal  of  optical  oceanographers. 
These  inherent  optical  properties1  (IOP’s)  affect  the 
ocean  surface  color,  the  transfer  of  heat  to  the  upper 
ocean,  the  transmission  of  photosynthetically  avail¬ 
able  radiation  through  the  water  column,  and  under¬ 
water  visibility.  The  value  of  a  is  also  used  in 
models  that  predict  phytoplankton  growth  rate  and 
ocean  primary  production,2  and  in  situ  measure¬ 
ments  of  a  and  bb  are  necessary  to  validate  remote 
sensing  algorithms  designed  to  monitor  IOP’s  on  a 
global  scale. 

A  common  method  for  determining  a  is  spectropho- 
tometric  analysis  of  discrete  water  samples.3  How¬ 
ever,  this  method  is  time-consuming,  has  a  limited 
sampling  rate,  and  is  subject  to  errors.  Alterna¬ 
tively,  a  can  be  determined  from  in  situ  natural  light 
measurements.  Most  simply,  a  can  be  determined 
from  simultaneous  in  situ  monochromatic  irradiance 
and  monochromatic  scalar  irradiance  measurements 
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with  the  Gershun  law4;  however,  monochromatic  sca¬ 
lar  irradiance  detectors  are  not  yet  readily  available. 
Instead,  estimates  of  a  have  been  made  from  near¬ 
surface  irradiance  measurements  in  conjunction  with 
measurements  of  remote  sensing  reflectance5  or  esti¬ 
mates  of  the  downward  mean  cosine  of  the  radiance 
distribution.6  However,  the  former  requires  addi¬ 
tional  above-surface  measurements  and  accurate  em¬ 
pirical  correlations,  whereas  the  latter  assumes  the 
downward  mean  cosine  does  not  change  significantly 
within  the  surface  layer,  and  both  methods  are  sus¬ 
ceptible  to  wave-induced  fluctuations  and  to  ship 
shadow. 

Reflecting-tube  instruments  make  it  possible  to  ob¬ 
tain  small-volume  in  situ  estimates  of  a  and  c  =  a  + 
b, 7  where  the  beam  attenuation  coefficient  c  is  the 
inverse  of  the  mean  free  path  of  a  photon.  Because 
reflecting-tube  instruments  are  subject  to  scattering 
errors,  they  use  a  small  sampling  volume  to  minimize 
these  errors,  and  as  a  result  these  instruments  can 
break  up  or  fail  to  collect  large  optically  active  aggre¬ 
gates  (marine  snow),  can  give  high  readings  that  are 
due  to  rare  events  of  large  particles  entering  the  sens¬ 
ing  area,  and  can  have  difficulty  detecting  low  con¬ 
stituent  concentrations. 

In  contrast,  IOP  estimation  from  natural  irradi¬ 
ances  can  be  obtained  from  large-volume  measure¬ 
ments,  and  therefore  small  concentrations  of 
constituents,  both  large  and  small,  can  be  detected.8 
Another  advantage  of  calculating  IOP’s  from  irradi¬ 
ance  measurements  is  that  it  enables  one  to  obtain 
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the  water  properties  and  light  field  from  the  same 
instrument.  At  the  very  least,  large-volume  mea¬ 
surements  can  be  correlated  with  small-volume  mea¬ 
surements  to  improve  the  confidence  in  these 
estimates.  Perhaps  the  primary  advantage  of  the 
large-volume  methods  is  that  much  irradiance  data 
have  already  been  collected  and  archived  that  can  be 
reanalyzed  with  new  algorithms  for  the  estimation  of 
a  and  bb. 

Here  we  evaluate  the  estimation  of  a  and  bb  from 
only  in  situ  profiles  of  the  upward  and  downward 
irradiances  Eu(z)  and  Ed(z)  at  geometric  depths  2. 
We  note  that  Gordon  and  Boynton9  demonstrated 
that  good  estimates  of  b  are  not  possible  from  only 
irradiance  measurements,  and  for  this  reason  we  fo¬ 
cus  on  estimating  a  and  bb.  Our  approach  is  to  de¬ 
termine  a  and  bh  through  the  determination  of  the 
reflectance  R{z )  and  the  downward  diffuse  attenua¬ 
tion  coefficient  Kd(z).  The  values  of  R(z)  and  Kd[z) 
are  used  to  estimate  the  IOP’s  R^  and  K,_,  which  are 
the  values  far  from  the  surface  of  R(z)  and  Kd(z), 
respectively.  Given  a  specific  scattering  phase  func¬ 
tion  [3,  there  are  unique  correlations  between  the  val¬ 
ues  of  and  K„  and  those  of  a  and  bh  that  can  be 
derived  from  the  radiative  transfer  equation. 

The  relevant  equations  of  radiative  transfer  are 
introduced  in  Section  2.  Estimation  of  R  y  and  K,  in 
deep,  homogeneous  waters  is  considered  in  Section  3, 
and  the  method  of  calculating  a  and  bh  from  K,  and 
R  ,_  is  presented  in  Sections  4  and  5.  The  importance 
of  selecting  an  appropriate  scattering  phase  function 
is  investigated  in  Section  6,  and  a  simplified  algo¬ 
rithm  that  is  independent  of  the  scattering  phase 
function  is  evaluated  in  Section  7.  In  Sections  8  and 
9  we  consider  cases  in  which  the  water  is  optically 
shallow  and  inhomogeneous,  respectively. 


tering  asymmetry  factor.  The  backscattering 
coefficient, 


bh  =  b 


P(fL  l)dp, 


(3) 


can  be  calculated  from10 


bh  =  bjb  =  (1/2) 


1  “  2  (2?l  +  1  )fn 


i  odd 


(4) 


once  the  fn  coefficients  are  specified.  The  integral 
factors  in  Eq.  (4)  can  be  calculated  numerically  from 
the  recursion  relationship  {n  +  1)  JV  P„(p)dp  = 
~(n  -  2)  Jo1  Pra_2(p)dp,  starting  with  V  Px(p)dp  = 
0.5. 

The  irradiance  reflectance  R(z)  is 


R(z)  =  Eu(z)/Ed(z),  (5) 

where  Eu(z)  and  Ed(z)  are  the  upward  and  downward 
irradiances: 


Eu(z)  =  |  |p-|£(z,  P-)dp, 

(6) 

Ed{z)  =  [  pL(z,  p)dp. 

(7) 

The  downward  irradiance  diffuse  attenuation  coeffi¬ 
cient  Kd  is  defined  by 


Kd(z)  = 


1  d Ed{z)  d  In [Ed{z)\ 


Ed(z)  dz 


d z 


(8) 


2.  Basic  Equations 

The  integrodifferential  transfer  equation  for  waters 
with  homogeneous  optical  properties  and  no  internal 
sources  is 


p j)L(z,  p )/dz  +  cL(z,  p)  =  b  J  (3(p,  p ')L(z,  p')dp', 

(1) 


Although  the  magnitudes  of  R{z)  and  Kjz)  near 
the  surface  depend  on  the  surface  illumination,  at 
large  depths  in  deep  homogeneous  waters  with  no 
internal  sources  the  values  of  R(z )  and  Kd(z)  ap¬ 
proach  asymptotic  values  P,_  and  K,,  respectively, 
that  are  IOP’s.  To  evaluate  R  ,,  and  Km  we  separate 
the  spatial  and  angular  dependencies  in  Eq.  (1)  with 
the  eigenmodes: 

L(z,  p)  =  c|4±v;,  pfexplTcz/r,)  (9) 


where  L(z,  pj  is  the  radiance  integrated  over  all  az¬ 
imuthal  directions  for  polar  angle  cos-1  p  with  re¬ 
spect  to  the  depth  z.  All  quantities  in  Eq.  (1)  are 
implicitly  a  function  of  wavelength.  The  azimuth- 
ally  integrated  scattering  phase  function  (3(p,  p')  is 
normalized  such  that  its  expansion  in  Legendre  poly¬ 
nomials  has  the  form 

1  M 

P(p,  P<' )  =  -  2  (2n  +  fo  =  1,  (2) 

^  n—0 

where  fn  are  the  expansion  coefficients,  P(1(p)  are  the 
Legendre  polynomials,  and  M  is  the  degree  of  scat¬ 
tering  anisotropy.  The  coefficient  fx  =  g  is  the  scat- 


and  use  Eq.  (2)  to  find  that  the  discrete  eigenfunc¬ 
tions  4> ( ±  v,- ,  p)  satisfy 


<K±v/,  p) 


2  (vj  +  p) 

M 

X  2  +  1  )fSn(±Vj)P,X[l),  Vj  >  1, 

n= 0 

(10) 


where  w0  is  the  single-scattering  albedo  «0  =  b/c. 
The  Chandrasekhar  polynomials11  gn  satisfy  the  re¬ 
cursion  formula 


ngjvj)  =  /i„_i vjg^ivj)  -  (n  -  1  )gn-2(vj),  GD 
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starting  with  g_±  =  0  and  g0  =  1,  where  lin  = 
(2 n  +  1)(1  —  oj(/„).  From  the  spherical  harmonics 
(P,v)  method12  with  N  odd  and  arbitrarily  large,  the 
positive  eigenvalues  Vj  are  approximately  the  roots  of 

gN+i&j)  =  0.  (12) 

With  this  formalism,  the  procedure  in  Appendix  A 
shows  that,  for  deep  homogeneous  waters,  R  ,_  and  Kx 
satisfy13 

=  Jo1  (|>(-vi,  pOp-dp, 

"  JV  4>(+  i’:.  p  )pdp  ' 

K„  =  c/v1,  (14) 

where  v1  is  the  largest  positive  eigenvalue.  Thus  Rx 
can  be  computed  directly  from  only  (3  and  oj0,  whereas 
K,_  can  be  determined  from  the  values  of  (3,  w0,  and  c. 
We  can  avoid  the  numerical  integration  in  Eq.  (13)  by 
utilizing  equations  for  computing  the  numerator  and 
denominator  of  Eq.  (13)  that  are  given  in  Ref.  14  [Eq. 
(22)]. 

3.  Estimation  of  and  K„  in  Deep,  Locally 
Homogeneous  Waters 

Given  the  upward  and  downward  irradiance  mea¬ 
surements  at  arbitrary  depths  2,  R(z)  is  calculated 
directly  from  Eq.  (5).  Values  of  Kd{z)  are  calculated 
from  finite  differences  of  ln[Ed(z)]  with  respect  to  z, 
an  approximation  that  is  exact  at  large  depths  where 
ln[Ed(z)]  varies  linearly  with  z.  Thus  the  calcula¬ 
tion  of  R(z)  at  a  specified  depth  requires  the  mea¬ 
surement  or  interpolation  of  Eu(z)  and  Ed(z)  at  that 
depth,  whereas  that  of  Kd(z)  requires  at  least  the 
measurement  of  Ed(z)  at  two  depths. 

In  optically  deep,  source-free,  and  homogeneous 
waters,  the  vertical  profiles  of  R(z)  and  Kd(z)  ap¬ 
proach  the  asymptotic  values  R  ,_  and  K,  that  we  seek. 
The  simplest  approach  to  estimate  R  ,  and  Kx  is  to 
calculate  R{z)  and  Kd[z),  preferably  at  large  depths, 
and  take  Rx  =  R(z)  and  K,  =  Kd{z).  We  refer  to  this 
as  the  asymptotic  method  (AM)  since  the  approxima¬ 
tion  is  exact  in  the  asymptotic  regime.  The  accuracy 
of  this  method  depends  on  the  degree  to  which  the 
light  field  differs  from  the  asymptotic  field  and  on  the 
noise  in  the  irradiance  measurements. 

If  measurement  noise  were  negligible,  the  AM 
would  be  essentially  exact  at  sufficiently  large 
depths.  For  illustration,  consider  the  case  in  which 
the  incident  radiation  is  composed  of  70%  direct  beam 
and  30%  diffuse  skylight  and  the  cosine  of  the  angle 
of  the  direct  beam  is  p0  =  0.866.  The  corresponding 
upper-boundary  condition  just  above  the  sea  surface 
can  be  modeled  as  a  superposition  of  direct  sunlight, 
defined  with  a  Dirac  delta  function,  and  diffuse  sky¬ 
light: 

L(0~,  p)  =  0.78(p  -  0.866)  +  0.3,  0  <  p  <  1.  (15) 

Let  co0  =  0.7  and  let  the  scattering  be  characterized 
by  the  Henyey-Greenstein  scattering  phase  func¬ 
tion15  (3HG,  for  which  fn  =  g",  with  the  scattering 


Fig.  1.  Diffuse  attenuation  coefficient  and  irradiance  reflectance 
profiles  from  simulated  irradiance  data  (o>0  =  0.7,  g  =  0.85).  The 
values  of  R  are  multiplied  by  10. 


asymmetry  factor  g  =  0.85.  We  take  the  200th- 
order  scattering  anisotropy  (M  =  200)  and  an  index  of 
refraction  of  1.34.  From  Eqs.  (13)  and  (14),  we  cal¬ 
culate  Rx  =  0.04119  and  KJc  =  0.4707.  The  value 
of  Kr_  is  normalized  by  c  so  that  the  result  depends 
only  on  w0  and  the  fn  values.  Let  the  vertical  dis¬ 
tance  be  measured  in  optical  depths  t  =  cz  and  take 
the  water  to  be  very  deep  (50  optical  depths)  with  a 
purely  absorbing  bottom.  As  shown  in  Fig.  1,  the 
profiles  of  R{ t)  and  KdJ)/c  converge  to  their  asymp¬ 
totic  values.  The  profiles  shown  in  Fig.  1  were  com¬ 
puted  with  the  discrete  ordinates  radiative  transfer 
code  disortb,  16,17  which  takes  into  account  the  index 
of  refraction  mismatch  at  the  surface.  The  magni¬ 
tudes  of  R( t)  and  Kd(i)/c  below  t  =  2  each  vary  over 
a  range  of  approximately  6%  and  below  t  =  4  vary  by 
less  than  3%.  This  small  range  of  R( t)  and  Kd( t) 
values  away  from  the  surface  is  typical  and  aids  in 
the  estimation  of  KJc  and  Rx. 

Table  1  shows  the  error  in  values  of  R  r_  and  Kf 
predicted  by  the  AM  from  the  simulated  data  shown 
in  Fig.  1.  Irradiance  values  at  1  optical  depth  ver¬ 
tical  spacing  were  used  to  calculate  Kd( t).  The  er¬ 
rors  in  the  predictions  decrease  monotonically  with 
depth.  However,  note  that  predictions  from  mea¬ 
surements  below  4  optical  depths  are  accurate  to 
within  a  few  percent.  This  indicates  that  even 


Table  1.  Predictions  of  Kx  and  Rx  from  the  Simulated  Kd{ t)  and  R{ t) 
Profiles  of  Fig.  1  versus  the  Depth  of  the  Deepest  Irradiance 
Measurement 


T 

Error  (%) 

AM“ 

EM" 

Rx 

KJc 

Rx 

KJc 

3 

4.0 

4.5 

0.78 

1.3 

4 

2.5 

2.8 

0.32 

0.48 

5 

1.6 

1.8 

0.12 

0.19 

7 

0.68 

0.75 

0.012 

0.018 

“Estimates  are  made  from  the  deepest  values  of  R( t)  and  Kd( t) 
(AM)  and  from  extrapolation  with  an  exponential  model  (EM)  for 
Kd( T)  and  R( t). 
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though  the  AM  is  only  theoretically  exact  at  large 
depths,  it  can  produce  good  estimates  from  relatively 
shallow  measurements. 

In  practice,  irradiance  measurements  contain 
noise.  Consider  random  noise  r  in  the  irradiance 
measurements  that  is  proportional  to  the  irradiance 
magnitude.  The  measured  reflectance  Rm{z)  is 

Rm(z )  =  I  1  +  y^W),  (16) 

where  A(z)  is  the  true  reflectance  and  ru  and  rd  are 
the  uncorrelated  random  noise  in  the  upward  and 
downward  irradiance  measurements,  respectively, 
after  any  smoothing  or  averaging  of  the  data.  The 
measured  downward  diffuse  attenuation  coefficient  is 

Kdm(z )  =  Kd(z)  -  d[ln(l  +  r„)]/dz.  (17) 

When  the  AM  is  used,  the  relative  error  of  R  ,.  is 
independent  of  R(z)  whereas  that  of  K„  is  propor¬ 
tional  to  \AzKJz)]  1,  where  A z  is  the  vertical  spacing 
between  Kd[z)  values.  For  example,  if  the  noise  rd 
and  ru  is  normally  distributed,  then  the  standard 
deviations  of  ( ru  -  rd)/(  1  +  rd)  and  [ln(l  +  rd  i+1)  - 
ln(l  +  rdi)\  are  both  approximately  1.4s  if  s  is  the 
standard  deviation  of  both  ru  and  rd.  In  such  a  case, 
the  relative  errors  in  A,.  and  K,  at  large  depths  are 
1.4s  and  l.As/\Kd(z)kz],  respectively. 

In  addition  to  the  AM  for  estimating  R  ,  and  K,  , 
one  can  use  analytical  approximations  to  R(z)  (Ref. 
18)  and  Kd(z)  (Ref.  19)  given  by 

R(z)  ~  A«,  +  [R(zr )  -  RJexp[-Sl>(z  -  zr)\  z  >  zr, 

(18) 

Kd(z)  =  K„  +  \Kd(zr)  -  K^expl-^iz  -  zr)\  z  >  zr, 

(19) 

where  zr  is  some  reference  depth  and  the  IOP  can 
be  calculated20  from  9?  =  c(v2  1  —  v  A '  ) ,  where  v2  is 
the  second  largest  positive  eigenvalue  from  Eq.  (12). 
Because  values  of  R(z)  and  Kd(z)  are  predicted  for 
depths  below  those  of  the  measurements,  we  refer  to 
this  method  as  the  extrapolation  method  (EM).  At 
large  depths  the  exponential  term  in  Eqs.  (18)  and 
(19)  becomes  negligible  and  EM  estimates  become 
equal  to  those  from  the  AM.  From  values  of  R(z )  at 
depths  z0,  zlf  and  z2,  with  zL  =  ( z0  +  z2)/ 2,  Rx  can  be 
obtained  from  Eq.  (18)  as20 

R(z0)R(z2)-R2(Zl) 

00  R(z0)+R(z2)-2R(z1)' 

An  analogous  equation19  holds  for  the  determination 
of  K„  from  Eq.  (19). 

Estimates  of  R  ,  and  K,  from  the  EM  for  the  noise- 
free  data  in  Fig.  1  are  given  in  Table  1.  Irradiance 
values  at  1  optical  depth  vertical  spacing  were  used, 
and  the  results  are  presented  as  a  function  of  the 
depth  of  the  deepest  irradiance  value  used.  All  pre¬ 
dictions  are  within  0.8%  of  the  true  value  and,  as 
expected,  improve  with  increasing  depth.  This  indi¬ 
cates  that  R  ,  and  K,  can  be  estimated  from  relatively 
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shallow  measurements,  at  least  when  there  is  no  sim¬ 
ulated  noise.  For  example,  at  t  =  3,  R  ,_  calculated 
from  the  EM  is  within  0.78%  of  the  true  value, 
whereas  A,_  estimated  by  the  AM  is  in  error  by  4%. 
Similarly,  at  t  =  4,  K,  calculated  from  the  EM  is 
within  0.48%  of  the  true  value,  whereas  K „  estimated 
by  the  AM  is  in  error  by  2.8%. 

When  noise  is  present,  however,  estimation  of  A  , 
and  K,  from  the  EM  can  give  improved  results  over 
the  AM  only  if  the  noise  in  R(z)  and  Kd(z )  is  much 
smaller  than  the  exponential  terms  in  Eqs.  (18)  and 
(19)  at  the  depths  of  the  measurements.  Thus  the 
extrapolation  method  may  be  superior  near  the  sur¬ 
face,  but  it  is  inferior  to  the  AM  at  depths  where  the 
true  R(z)  and  Kd{z)  vary  only  slowly  with  depth.  Be¬ 
cause  our  emphasis  is  on  deeper-water  applications, 
the  EM  is  not  examined  further. 

4.  Estimation  of  Fundamental  lOP’s  from  R „  and  K „ 

The  parameters  a,  b,  and  bb  are  fundamental  IOP’s. 
In  addition,  Gordon  et  al.21  showed  that  R( t)  at  the 
surface  is  nearly  proportional  to  the  ratio  bj (a  +  bb), 
and  many  remote  sensing  algorithms  are  designed  to 
measure  either  bb/a  or  bb/(a  +  bb).  We  choose  to 
define  the  ratio  G  =  bb/(a  +  bb)  as  the  (dimension¬ 
less)  Gordon  parameter  in  view  of  his  many  contri¬ 
butions  to  the  field  of  ocean  optics  and  his  observation 
of  the  importance  of  G  in  remote  sensing  applications. 

If  the  scattering  phase  function  (3  is  assumed,  then 
a  and  bb  can  be  determined  from  measurements  of  A, 
and  Kr  .  The  first  step  in  the  solution  of  this  inverse 
problem  is  to  calculate  w0  from  R  ,  and  (3,  or  from  A  ,_ 
and  g  if  the  phase  function  (3HG  is  assumed.  For 
example,  Fig.  2  shows  the  interdependence  of  A  /:,  co0, 
and  g  for  the  phase  function  fiHt;.  Although  in  a 
typical  forward  calculation  R  ,  depends  on  w0  and  g, 
in  this  inverse  problem  R,_  is  a  measured  quantity. 
The  value  of  oj0(Ax,  g)  can  be  found  from  an  iterative 
solution  of  Eq.  (13)  with  the  help  of  Eqs.  (10)— ( 12), 
and  in  the  process  the  value  of  v1  is  calculated  from 
Eq.  (12).  Next,  c  is  calculated  from  v1  and  the  mea¬ 
sured  value  of  K,  from  Eq.  ( 14),  a  and  b  are  calculated 
from  co0  and  c,  and  bh  is  computed  from  Eq.  (4)  with 
use  of  fn  =  g"  with  the  assumed  g.  Finally,  ratios 
such  as  bb/a  and  bb/(a  +  bb)  can  be  formed  from  a 
and  bb. 


Table  2.  Percent  Errors  in  the  Estimates  of  IOP’s  for  Selected  Values 


of  Percent  Errors  in  Rx 

and  Kx 

Error  (%) 

It, 

K, 

a 

b  and  bh 

bb/a 

G 

5 

0 

-1.04 

3.32 

4.41 

4.06 

-5 

0 

1.08 

-3.38 

-4.41 

-4.08 

5 

5 

3.91 

8.49 

4.41 

4.06 

-5 

5 

6.14 

1.45 

-4.41 

-4.08 

5 

-5 

-5.99 

-1.84 

4.41 

4.06 

-5 

-5 

-3.97 

-8.21 

-4.41 

-4.08 

0 

5 

5.00 

5.00 

0 

0 

Table  2  gives  example  calculations  of  the  errors  in 
the  estimated  values  of  a,  bb,  bb/a,  and  G  obtained 
from  various  values  of  R  ,  and  Kr  when  the  true  IOP’s 
are  to0  =  0.7  and#  =  0.85.  These  computations  were 
performed  with  a  =  0.15  and  b  =  0.35  and  the  correct 
phase  function  and  by  an  iterative  search  for  the 
optimal  value  of  w0.  We  achieved  the  search  by  min¬ 
imizing  the  error  in  the  calculated  value  of  R  ,_  with 
respect  to  w0  using  Brent’s  method,22  which  is  a  com¬ 
bination  of  an  inverse  parabolic  interpolation  and  a 
golden-section  search.  It  can  be  seen  that  the  errors 
in  the  calculations  of  the  fundamental  IOP’s  are  of 
the  same  order  of  magnitude  as  the  errors  in  Rrj0  and 
K,  .  Estimates  of  a  are  best  when  the  errors  in 
and  K,  have  the  same  sign,  whereas  estimates  of  bh 
are  best  when  the  errors  in  R,  and  K„  have  the 
opposite  sign.  As  can  be  seen  from  Table  2,  bb/a  and 
G  do  not  depend  on  K„  since  bb/a  =  56oo0(l  -  w0)_1 
and  G  =  [1  +  (1  —  co0)/(co0b6)]“1  and  w0  is  indepen¬ 
dent  of  K 

In  practice,  fi  is  not  well  known,  and  therefore  es¬ 
timates  of  co0  and  c  (and  of  a  and  bb)  will  contain 
errors  that  are  due  both  to  measurement  errors  in  A  , 
and  K,  and  to  the  error  in  the  assumed  \i.  Table  3 
shows  the  percent  errors  in  estimates  of  a  and  bb 
calculated  with  the  same  iterative  solution  code  as  for 
Table  2  but  with  unknown  g  and  for  R  ,  and  K,  that 
were  computed  for  the  indicated  values  of  co0  and  g. 
For  the  assumed  g  =  0.85,  the  values  of  a  and  bh  are 
underestimated  when  the  true  g  <  0.85,  are  exact 
when g  =  0.85,  and  are  overestimated  when g  >  0.85. 
For  this  large  range  of  g,  the  worst  estimate  of  a  was 


Table  3.  Percent  Errors  in  the  Estimates  of  a  and  bb  for  Waters  with 
Given  Values  of  g  and  co0“ 


g 

Error  (%) 

“o  = 

0.5 

“o  - 

0.7 

O>0  = 

0.9 

a 

bb 

a 

bb 

a 

bb 

0.75 

-1.6 

-8.9 

-1.2 

-9.6 

-0.52 

-9.8 

0.80 

-0.83 

-4.4 

-0.65 

-4.9 

-0.30 

-5.1 

0.90 

0.86 

3.7 

0.78 

4.6 

0.45 

5.3 

0.95 

1.6 

6.1 

1.7 

7.9 

1.2 

10 

“The  estimates  were  obtained  assuming#  =  0.85  and  using  the 
and  Kx  values  for  the  indicated  g  and  u>0. 


Fig.  3.  Normalized  sensitivity  coefficient  of  u>0  with  respect  to  Rx 
for  g  of  Phg- 


in  error  by  only  1.5%,  whereas  those  for  bb  were 
generally  within  10%. 

To  assess  more  thoroughly  the  accuracy  of  predict¬ 
ing  fundamental  IOP’s  from  R  ,_  and  K,,  in  Section  5 
we  examine  the  sensitivity  coefficients  that  quantify 
the  extent  to  which  errors  in  R  r  i  K,_,  and  g  affect  the 
estimates  of  the  fundamental  IOP’s. 

5.  Sensitivity  Coefficients 

Normalized  sensitivity  coefficients  express  the  ratio 
of  the  relative  error  in  an  output  (e.g.,  a)  to  a  small 
relative  error  in  an  input  (e.g.,  Ra 0).  In  the  develop¬ 
ment  of  equations  for  sensitivity  coefficients  in  this 
section,  it  is  assumed  that  the  phase  function  can  be 
expressed  as  a  function  of  a  single  parameter,  such  as 
the  scattering  asymmetry  factor  g  of  the  Henyey- 
Greenstein  phase  function  that  is  used  for  the  nu¬ 
merical  calculations. 

A.  Sensitivity  Coefficients  for  oo0  and  c 

From  the  iterative  search  method  discussed  in  Section 
4,  one  can  numerically  solve  for  (R.,J w0)(dw0/ dR  j  and 
(g/(o0)(d(x)0/dg).  These  normalized  sensitivity  coeffi¬ 
cients  quantitatively  express  the  sensitivity  of  the  es¬ 
timates  of  co0  to  errors  in  the  measurement  of  R  ,_  and 
in  the  guess  for  g,  respectively.  Furthermore,  these 
coefficients  are  required  to  calculate  the  sensitivity 
coefficients  for  the  other  fundamental  IOP’s.  As 
shown  in  Appendix  B,  these  sensitivity  coefficients  can 
also  be  expressed  in  terms  of  sensitivity  coefficients  for 
the  forward  problem.  Because  the  forward  problem  is 
much  easier  to  compute  than  the  inverse  problem,  iR.,J 
w0)(da)0 /dRoo)  and  (g/ oj0)(V')oj0/ i)g)  were  computed  with 
Eqs.  (Bl)  and  (B2). 

In  Fig.  3  (R,yJ gj0)  ((5to0/ <)R,y)  is  shown  as  a  function  of 
co0.  For  0.75  <  g  <  0.95,  this  normalized  sensitivity 
coefficient  varies  roughly  linearly  from  approxi¬ 
mately  0.7  for  w0  =  0.2  to  0.05  for  w0  =  0.99.  This 
indicates  that  estimates  of  w0  from  R  ,  are  moderately 
insensitive  to  small  errors  in  R  ,  and  are  least  sensi¬ 
tive  where  the  absorption  relative  to  the  beam  atten¬ 
uation  is  lowest.  The  sensitivity  coefficient  itself  is 
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Fig.  4.  Normalized  sensitivity  coefficient  of  oo0  with  respect  to  g  of 
Phg- 


Fig.  6.  Normalized  sensitivity  coefficient  of  c  with  respect  to  g  of 
Phg- 


only  weakly  dependent  on  g  but  is  highest  for  small  g 
when  co0  >  0.5  and  for  large  g  when  co0  is  small. 
Unfortunately,  as  can  be  seen  in  Fig.  4,  estimates  of 
w0  are  much  more  sensitive  to  g  than  to  Rx.  The 
sensitivity  of  co0  to  g  is  highly  dependent  on  both  g 
and  oj0,  and  it  is  highest  for  large  g  and  for  small  co0. 
Thus  the  sensitivities  of  w0  to  errors  in  both#  and  Rx 
are  lowest,  and  therefore  estimates  of  w0  will  be  best, 
when  w0  is  high.  For  w0  =  0.7  and  g  =  0.85,  for 
example,  a  10%  uncertainty  in  g  can  result  in  approx¬ 
imately  an  18%  uncertainty  in  o»0,  which  would  be 
unacceptably  large.  Much  greater  errors  may  result 
when  g  is  very  high  (g  >  0.9),  which  is  typical  for  the 
Petzold  phase  functions.23  Because  we  wish  to  im¬ 
plement  this  method  of  estimating  IOP’s  without 
knowledge  of  g,  the  high  sensitivity  of  co0  to  g  places 
a  serious  limitation  on  our  ability  to  estimate  co0. 
However,  some  other  IOP’s  calculated  from  oj0,  in 
particular  a,  are  far  less  sensitive  to  g. 

In  the  inverse  problem,  the  largest  eigenvalue  v±  is 
a  function  of  oj0(i?.,_,  g)  and  g  and  therefore  can  be 
written  alternatively  as  a  function  of  the  independent 


variables  R,  and  g.  The  beam  attenuation  coeffi¬ 
cient  is  found  from 


c  =  vt(Rx,  g)K,,  (21) 


but  henceforth  the  dependence  of  iq  and  co0  on  R  ,_  and 
g  will  not  be  denoted  explicitly.  Note  that  K,  is  a 
measured  and  therefore  independent  quantity.  One 
can  evaluate  the  normalized  sensitivity  coefficients 
for  c  from  Eq.  (21): 


K  ,  dc 


(22) 


Rx  dc  Rx  dv1  Rx  dv1 
c  dRx  c  ”  dRx  v1  dR„  ’ 


g  dc  =  gKx  dig  =  g  dig 
eg  c  dg  vi  dg  ' 


(24) 


Fig.  5.  Normalized  sensitivity  coefficient  of  e  with  respect  to  Rx 
for  g  of  Phg- 


These  are  expressed  in  terms  of  normalized  sensitiv¬ 
ity  coefficients  (R.Jv1)(dv1/dRx)  and  (g/v1)(dv1/dg), 
which  must  be  computed  either  from  an  iterative 
solution  method  or,  more  easily,  from  Eqs.  (B3)  and 
(B4).  The  magnitude  of  (R,Jc)(dc/dRg>  shown  in 
Fig.  5,  is  typically  less  than  unity,  except  for  high  w0 
and  low  g,  and  is  especially  low  when  w0  is  low. 
However,  the  sensitivity  of  c  to  g,  shown  in  Fig.  6,  is 
similar  in  magnitude  to  that  of  w0  to  g,  except  that  it 
is  lowest  for  low  values  of  co0  and  highest  for  high 
values  of  w0. 

B.  Sensitivity  Coefficients  for  a,  bb,  bja,  and  G 
Because 

a  =  c(l  —  w0)>  da/dwo  =  — c,  da/dc  =  (1  —  w0)»  (25) 

b  =  cto0,  db/du0  =  c,  db/dc  =  w0,  (26) 
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the  normalized  sensitivity  coefficients  for  a  can  be 
determined  directly  from  those  for  w0  (Figs.  3  and  4) 
and  c  [Eqs.  (22)— (24)]: 


K,  da  K,  da  dc 
a  dK,  a  dc  dK, 


—  (1  -  w0)v,  =  1, 
a 


(27) 


R 

da 

R, 

da  c 

R-,  dc  \ 

da 

Wo 

(R, 

dw0\ 

— 

=  - 

— 

— 

+  — 

— 

a 

dR, 

a 

dc  R „ 

\  C  dR., 

O 

3 

<■£> 

R, 

\W0 

dR„)\ 

R,  dc 
c  dR. 


w0 


Wo, 


du>0 
w0  dR.,. 


(28) 


gda 
a  dg 


g  dac^lgdc\+da^u>olg_du>o 
a  [dc  g  \c  dg)  dco0  g  \w0  dg 
fgdc]  _  /  w0  \/g  dw0\ 

\C  \1  —  w0/ \w0  dg  ) 


(29) 


Similarly,  the  normalized  sensitivity  coefficients  for  b 
are 


K,  db 

R.  db 

b  dR  ., 


K, 

—  WqV!  =  1, 


(30) 


R,  db  c  R„  dc  \  db  w0 
b  dc  R „  \  c  dR.,)  dw0  R, 
(Rx  dc  \  (R-,  dw0\ 

\c  dRj  +  \u0dRj’ 


R,  dw0\ 
,w0  dR„) 


(31) 


gdb^_g  db  C  Ig  dc\  |  db  (jQq  /  g  dfa)Q 
b  dg  b  [dc  g  \c  dg)  dco0  g  \w0  dg 

=  (g  +  /  g  ^Op\ 

\C  dg)  \CD0  dg  )  ‘ 


(32) 


Thus,  all  else  being  equal,  the  percent  error  in  a  and 
b  that  is  due  to  an  error  in  K,  is  equal  to  the  percent 
error  in  K,,  The  normalized  sensitivity  coefficients 
for  b  to  R, o  and  g  are  equal  to  the  sum  of  those  for  c 
and  w0  to  R,  and g,  respectively,  whereas  the  expres¬ 
sions  for  the  normalized  sensitivity  coefficients  for  a 
to  R  ,  and  g  are  similar  to  those  for  b  except  that  for 
a  the  coefficients  involving  co0  are  scaled  by  —  co0/ (1  - 
w0)  =  -b/a.  Note  that  because  dc/ dR„  and  <ioj0/ dR  , 
have  the  same  sign,  the  sensitivity  coefficients  for  a 
are  less  than  those  for  b,  which  is  consistent  with  the 
observation  in  Section  4. 

The  coefficient  (g/a)(da/dg)  is  shown  in  Fig.  7. 
The  absolute  value  of  the  magnitude  of  the  coefficient 
is  very  low  (<0. 16)  for  co0  <  0.99.  It  is  highly  depen¬ 
dent  on  w0,  with  the  largest  magnitude  corresponding 
to  moderate  values  of  «0.  The  coefficient  (R  ,/a)(da/ 
dR,)  is  shown  in  Fig.  8.  This  coefficient  also  is  rel¬ 
atively  small,  with  typical  values  ranging  from  —0.1 
to  —0.6  for  oj0  <  0.9.  Therefore,  even  though  a  is 
calculated  from  co0  and  c,  each  of  which  are  poorly 
estimated  from  R,  and  K„  a  can  be  calculated  quite 
well.  As  expected,  the  sensitivity  of  b  to  g  as  calcu¬ 
lated  from  Eq.  (32)  is  very  large,  especially  for  large 


Fig.  7.  Normalized  sensitivity  coefficient  of  a  with  respect  to  g  of 
Phg- 


values  of  g,  and  therefore  reasonable  estimates  of  b 
when  p  is  not  known  are  not  possible. 

The  normalized  sensitivity  coefficient  of  bb  to  g  is 
the  sum  of  those  of  b  to  g  and  bh  to  g: 


g  dbb  _gdb  ^  g  dbb 
bb  dg  b  dg  bh  dg  ' 


This  introduces  a  fifth  sensitivity  coefficient  {g/ 
bb)(dbb/dg)  that  must  be  computed  numerically.  A 
plot  of  (gdbb)/{bhdg)  is  shown  in  Fig.  9.  The  terms 
(g/b)(db/dg)  and  (g/bb)(dbh/dg)  tend  to  cancel  out  be¬ 
cause  the  signs  are  different,  so  estimates  of  bb  are 
moderately  insensitive  to  g  [0.4  <  \(gdbb)/(bbdg)\  < 
0.9],  indicating  that  estimates  of  bb  will  be  reasonable 
even  if  g  is  not  well  known.  Because  bh  =  bbb,  the 
normalized  sensitivity  coefficient  of  bb  to  K,  is  unity, 
whereas  that  of  bb  to  R  ,  is  identical  to  that  for  b  given 
in  Eq.  (31)  and  is  shown  in  Fig.  10. 


Fig.  8.  Normalized  sensitivity  coefficient  of  a  with  respect  to  R„ 
for  g  of  (1HG. 
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Fig.  9.  Normalized  sensitivity  coefficient  of  bh  with  respect  to  g  of 
Phg- 


Fig.  11.  Normalized  sensitivity  coefficient  of  G  with  respect  to  Rx 
for  g  of  (lHG. 


Recall  that  both  bb/a  and  G  are  independent  of  KA. 
The  normalized  sensitivity  coefficient  of  bb/a  to  g  is 


ficients  of  b,J  a  and  G  to  R,rj  are  in  the  same  form  as 
those  to  g: 


g  d(bb/ a) 

(bb/a)  dg 


g  dbb  g  da 

bh  dg  a  dg 

g  dbb  |  1  g  d(0 o 

bh  dg  (1  -  co0)  w0  dg  ’ 


(34) 


and  the  sensitivity  coefficient  of  G  to  g  is  proportional 
to  that  of  ( bb/a )  to  g: 


Ro=  d(bb/a)  R *  dbb  1  R™  dw0 
( bb/a )  dR„  bb  dR<»  (1  —  w0)  co0  dR„  ’ 


(36) 


R„  dG  a 
G  dR„  a  +  bb 

The  sensitivities  of  G  are  shown  in  Figs.  11  and  12. 


R*  d(bb/a) 
(■ bb/a )  dRoo 


(37) 


g  dG  a 
G  dg  a  +  bb 

Because  dbb/dg  and  da/dg  have  the  same  sign  and 
\(gdbb)/(bbdg)\  >  \(gda)/(adg)\,  then  from  Eq.  (34)  the 
magnitude  of  the  sensitivity  of  bb/a  to  g  is  less  than 
that  of  bb,  and  because  a/ (a  +  bb)  <  1,  then  from  Eq. 
(35)  G  is  even  less  sensitive  to  g  than  is  bb/a,  espe¬ 
cially  for  large  co0.  The  normalized  sensitivity  coef- 


6.  Choice  of  the  Phase  Function  Model 

Equations  (29)  and  (32)  hold  for  any  model  of  the 
phase  function  (3  that  can  be  specified  by  a  single 
parameter  g.  In  Figs.  7  and  9  the  one-term  Henyey- 
Greenstein  function  was  used  because  of  its  simplic¬ 
ity.  Other  models  could  be  used,  however.  For 
example,  phase  functions  obtained  by  a  combination 
of  the  pure  seawater  phase  function  and  the  Petzold 
particle  phase  function24  can  also  be  specified  by  the 


g  d(bb/a) 


(bh/a)  dg 


(35) 


Fig.  10.  Normalized  sensitivity  coefficient  of  b  or  bb  with  respect  Fig.  12.  Normalized  sensitivity  coefficient  ofG  with  respect  tog  of 
to  R.,_  for  g  of  Phg-  (W 
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value  of  g  or,  more  commonly,  the  chlorophyll  con¬ 
centration. 

Even  if  the  single  parameter  defining  the  phase 
function  model  is  known,  errors  may  be  introduced 
into  the  estimates  of  a  and  bb  if  the  shape  of  the  true 
phase  function  differs  from  that  predicted  by  the 
model.  For  example,  Fig.  13  shows  four  phase  func¬ 
tions  that  are  all  characterized  by  g  =  0.85  but  that 
have  very  different  shapes  in  the  backscattering  di¬ 
rections.  These  were  generated  from  the  two-term 
Henyey-Greenstein  model,25  which  is  a  linear  com¬ 
bination  of  two  one-term  Henyey-Greenstein  func¬ 
tions: 


P(gi,  g2,  «)  =  «P(gi)  +  (1  -  a)P(g2)>  (38) 

with  parameters  gv  g2,  and  a,  for  0  <  a  <  1.  Selection 
of  a  negative  g2,  for  example,  yields  an  enhanced  back- 
scattering  effect.  Shown  are  a,  the  one-term  function 
( g1  =  0.85,  a  =  0);  b,  a  monotonic  function  with  more 
backscattering  (gq  =  0.88,  g2  =  -0.062,  a  =  0.97);  c, 
the  function  whose  first  three  moments  match  those  of 
the  water-Petzold  particle  mixture  model24  (gt  =  0.90, 
g2  =  -0.26,  a  =  0.96);  and  d,  a  case  with  extremely 
high  backscattering  ( g1  =  0.90,  g2  =  -0.64,  a  =  0.97). 

The  phase  functions  in  Fig.  13  were  used  to  deter¬ 
mine  the  effect  of  the  amount  of  backscattering  in  the 
assumed  phase  function  on  the  estimates  of  the 
IOP’s.  Values  of  R  ,  and  K.,  were  first  calculated 
from  the  forward  problem  for  each  of  the  four  phase 
functions.  Then  the  values  of  a  and  bb  were  esti¬ 
mated  from  the  Rrj  and  K,  with  the  approach  in 
Section  4  assuming  either  phase  function  a  or  c.  Ta¬ 
ble  4  shows  the  percent  errors  in  the  estimated  IOP’s. 
Estimates  of  bb  were  several  times  more  sensitive  to 
the  assumed  (3  than  were  those  of  a.  Because  phase 
functions  a  and  b  are  decreasing  monotonically  in  the 
backscattering  direction,  estimates  of  bb  assuming 
phase  function  a  were  in  error  by  as  much  as  30% 
when  the  true  (3  exhibits  the  enhanced  backscatter- 


Fig.  13.  Two-term  Henyey-Greenstein  phase  function  model  for 
g  =  0.85  and  a,  g1  =  0.85,  a  =  0;  b,  g1  =  0.88,  g2  =  -0.062,  a  = 
0.97;  c,  gx  =  0.90,  g2  =  -0.26,  a  =  0.96;  and  d,  gx  =  0.90,  g2  = 
-0.64,  a  =  0.97. 


Table  4.  Percent  Errors  in  the  Estimates  of  IOP’s  for  io0  —  0.7  and  the 
Two-Term  Henyey-Greenstein  Phase  Functions  Shown  in  Fig.  13“ 


True  (3 

Error  (%) 

Assumed  |3  Phase 
Function  a 

Assumed  (3  Phase 
Function  c 

a 

K 

G 

a 

bb 

G 

a 

0.0 

0.0 

0.0 

5.8 

-52 

-49 

b 

-2.5 

-17 

-14 

3.1 

14 

9.5 

c 

-5.3 

-28 

-22 

0.0 

0.0 

0.0 

d 

-6.9 

-26 

-19 

-1.7 

2.8 

4.0 

“The  assumed  (3  is  either  the  one-term  phase  function  a  (a  =  0 
andgj  =  0.85)  or  the  two-term  function  c  =  0.90,  g2  =  -0.26, 
a  =  0.96).  The  true  phase  functions  used  in  the  forward  calcula¬ 
tion  of  and  K.,_  are  functions  a,  c,  b  (g1  =  0.88,  g2  =  -0.062,  a 
=  0.97),  and  d  (g,  =  0.90,  g2  =  -0.64,  a  =  0.97). 


ing  of  phase  functions  c  and  d.  Similarly,  assuming 
phase  function  c  when  the  true  phase  function  was  a 
or  b  led  to  large  errors.  Conversely,  when  phase 
function  c  was  assumed  and  the  true  (3  was  phase 
function  d,  which  also  exhibits  enhanced  backscatter¬ 
ing,  the  estimates  of  a  and  bb  were  in  error  by  only 
1.7%  and  2.8%,  respectively,  even  though  the  magni¬ 
tude  of  phase  function  d  at  180  deg  is  several  times 
larger  than  that  of  c. 

These  results  show  that,  to  obtain  good  estimates  of 
bb,  it  is  not  sufficient  for  the  assumed  phase  function 
to  be  characterized  by  an  appropriate  value  of  g;  the 
amount  of  backscattering  is  also  important.  In  par¬ 
ticular,  the  one-term  Henyey-Greenstein  function 
model  should  not  be  used  in  the  estimation  of  bh  in 
natural  waters;  the  Petzold  phase  function  or  a  two- 
term  Henyey-Greenstein  function  similar  to  c  would 
be  better. 


7.  Phase-Function-Independent  Algorithm 

Haltrin26  recognized  that  a  simple  model  for  an  ap¬ 
proximate  phase  function  |3(p,,  p/)  =  2 bh  +  2(1  — 
2&6)8(p  —  p')  conveys  much  of  the  information 
needed  to  describe  highly  forward  scattering.  For 
this  (3  he  derived  equations  for  R  ,_  and  KJc  from  a 
carefully  derived  two-flux  theory: 


a  _  (1  -  V(RJ2(1  +  4  +  R,J 

bb  4 Ry 


(39) 


KJc  =  (1  -  co0) 


1  +  G  ] 

1  +  2G  -  [G(4  +  5G)]1/2j 


1/2 


(40) 


Equation  (39)  can  be  used  to  estimate  a/bb  from  R„, 
and  G  can  be  calculated  from  G  =  (1  +  a/b 6)_1. 
Then  a  can  be  estimated  from  a  rearrangement  of  Eq. 
(40): 


1  +  G  ]"1/2 

1  +  2G  -  [G(4  +  5G)]1/2j  ’ 


(41) 


which  enables  one  to  also  obtain  bb.  The  primary 
difference  between  Haltrin’s  model  and  that  proposed 
in  Section  4  is  that  in  Haltrin’s  model  only  two  closed- 
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Table  5.  Percent  Errors  in  the  Estimates  of  a  and  bb  Obtained  from 
Haltrin’s  Model  for  c  =  1  and  the  One-Term  Henyey-Greenstein  Phase 
Function" 


g 

Error  (%) 

O)0 

=  0.5 

co0 

=  0.7 

O)0 

=  0.9 

a 

bb 

a 

bb 

a 

bb 

0.75 

2.2 

50 

2.5 

67 

1.4 

82 

0.80 

2.8 

51 

3.2 

71 

2.0 

90 

0.85 

3.2 

50 

4.0 

73 

2.7 

98 

0.90 

3.4 

44 

4.6 

69 

3.8 

100 

0.95 

2.8 

32 

4.5 

53 

5.2 

97 

“Values  of  and  Ka 0  were  calculated  for  the  given  a>0  andg,  and 

were  used  to  estimate  a  and  bb  from  Eqs.  (39)  and  (41). 


form  analytical  equations  are  needed  and  no  assump¬ 
tion  about  the  phase  function  is  required,  whereas  in 
the  method  proposed  in  Section  4  an  assumed  phase 
function  must  be  incorporated  into  the  iterative  so¬ 
lution  of  Eqs.  (10)-(13).  Although  Eqs.  (39)  and  (41) 
have  the  significant  advantage  in  that  they  are  easier 
to  implement  in  practice,  they  are  less  flexible  since 
they  do  not  admit  a  priori  information  about  the 
phase  function. 

Table  5  shows  the  percent  errors  in  estimates  of  a 
and  bb  obtained  from  Eqs.  (39)  and  (41)  for  the 
one-term  Henyey-Greenstein  phase  function  with 
0.75  <  g  <  0.95  and  0.5  <  co0  <  0.9.  As  for  Table 
3,  the  values  of  R  ,_  and  K,r  were  computed  for  the 
given  values  of  co0  and  g.  The  calculations  were 
performed  for  c  =  1;  however,  the  percent  errors  in 
a  and  bb  were  found  to  be  insensitive  to  the  value  of 
c.  The  errors  in  the  estimates  of  a  range  from  1.4% 
to  5.2%,  whereas  those  in  the  estimates  of  bb  range 
from  32%  to  100%.  The  errors  in  a  are  roughly 
twice  those  in  Table  3  obtained  with  the  approach 
described  in  Section  4,  whereas  the  errors  in  bb  are 
many  times  larger.  Similar  results  were  obtained 
with  the  San  Diego  Petzold  water  phase  function  for 
a  case  of  10-mg/m3  chlorophyll  concentration  and 
685-nm  light;  for  «0  =  0.5,  0.7,  and  0.9  the  percent 
errors  in  a  were  1.8,  2.6,  and  1.9,  respectively,  and 
the  percent  errors  in  bb  were  44,  61,  and  81,  respec¬ 
tively. 

The  large  errors  in  bb  should  be  expected  because  of 
the  high  sensitivity  of  bh  to  the  shape  of  the  backscat- 
tering  portion  of  (3  demonstrated  in  Section  6.  The 
degree  to  which  the  (3-dependent  approach  of  Section 
4  outperforms  Haltrin’s  approach  depends  on  how 
well  the  assumed  (3  matches  the  actual  scattering 
phase  function  of  the  water. 

8.  IOP  Estimation  in  Shallow  Waters 

In  shallow  water  where  the  entire  water  column  is 
in  the  euphotic  zone,  profiles  Riz)  and  Kd(z)  are 
affected  by  the  interaction  of  light  with  the  bottom. 
Although  the  effects  of  the  bottom  usually  do  not 
reach  as  far  into  the  water  column  as  do  surface 
conditions,  bottom  effects  can  cause  errors  in  IOP 
estimation.  If  the  water  is  so  shallow  that  signif¬ 


Fig.  14.  Estimation  of  Rx  in  water  of  5  optical  depths  with  a 
purely  absorbing  bottom  (Rb  =  0).  Shown  are  Rm  from  Eq.  (13), 
the  local  irradiance  ratio  R( t)  that  forms  the  asymptotic  model 
(AM),  and  the  depth-dependent  estimates  of  R,  from  the  deep- 
measurement  reflectance  model  (DMRM)  and  the  shallow- 
measurement  reflectance  model  (SMRM). 


icant  surface  and  bottom  effects  overlap,  then  Riz) 
and  Kd(z)  never  reach  R x  and  K,r,  respectively. 
However,  estimates  of  R,,_  and  can  still  be  made 
in  water  of  at  least  a  few  optical  depths.  In  this 
case,  Riz)  and  Kd(z)  near  the  surface  tend  toward 
R  ,  and  K,  with  increasing  depth  and  near  the  bot¬ 
tom  deviate  away  from  R  ,_  and  if  ,.  The  asymptotic 
method  for  estimating  R„,  for  example,  can  be  ap¬ 
plied  by  taking  R „  =  R(zm)  at  a  mid-water-column 
depth  zm  where  R(z)  is  a  maximum  or  minimum  or 
is  at  an  inflection  point. 

Figure  14  shows  R( t)  for  water  of  5  optical  depths 
with  a  purely  absorbing  bottom,  i.e.,  the  bottom  al¬ 
bedo  Rb  vanishes.  This  water  is  characterized  by  co0 
=  0.70,  g  =  0.85  and  fty  =  0.0412  based  only  on  the 
IOP’s.  From  top  to  bottom,  R( t)  increases  from  its 
surface  value  toward  the  value  of  R  ,_  and  then  de¬ 
creases  to  the  bottom  value  oiRb  =  0.  The  best  place 
to  apply  the  AM  in  this  case  is  at  the  maximum  R(t). 
Such  an  estimate  will  always  be  shy  of  the  true  value, 
but  it  is  difficult  to  know  by  how  much,  especially 
from  noisy  measurements.  Figure  15  shows  the 
same  water  as  in  Fig.  14,  but  with  a  Lambertian 
reflecting  bottom  of  Rb  =  0.2.  In  this  case  R{ t)  has 
an  inflection  point  near  the  depth  where  R(t)  =  Rm 
although  it  is  barely  detectable  because  of  the  large 
bottom  effects.  Because  Rb  :»  f?,„  values  of  R( t) 
deep  in  the  water  column  are  far  from  7?,. 

To  improve  our  ability  to  treat  shallow  waters,  two 
alternative  methods  for  calculating  R,  from  Ed(z) 
and  Eu{z)  near  the  bottom  are  derived  in  Appendix  C: 
the  deep-measurement  reflectance  model  (DMRM) 
and  the  shallow-measurement  reflectance  model 
(SMRM).  For  the  DMRM  one  must  make  measure¬ 
ments  of  Ed(z )  and  Eu(z)  at  two  depths,  subtract  and 
add  the  irradiances  at  each  depth,  square  the  results, 
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Fig.  15.  Estimation  of  R „  as  in  Fig.  14  for  a  Lambertian  bottom 
reflectance  of  Rb  =  0.2. 

and  evaluate  the  differences  at  the  two  depths.  The 
estimated  value  of  R  r_  follows  from 

/l"iU2  =  [EJz)-Eu{z)]\j 

\1  +Rj  [Ed(z)  +  £„(*)]  V  fi‘' 

Note  that  this  technique  for  estimating  Ry  can  be  es¬ 
pecially  useful  for  2  near  the  bottom  because  it  incor¬ 
porates  the  growing  eigenmode  g?1(:i:v1)exp(cz/v1)  (see 
Appendix  C)  that  is  needed  to  account  for  the  asymp¬ 
totic  contribution  of  the  bottom  boundary  condition. 
An  advantage  of  the  DMRM  is  that  the  bottom  albedo 
Rh  need  not  be  known,  but  a  disadvantage  is  that 
differences  of  noisy  irradiance  measurements  poten¬ 
tially  can  lead  to  large  errors. 

For  the  SMRM  one  needs  to  know  the  bottom  al¬ 
bedo  Rb  and  the  depth  of  the  water  column  zb.  From 
measurements  of  R(z)  at  arbitrary  geometric  depth  z, 
R  ,_  can  be  predicted  from 

R«,  =  {R(z)  -  Rb  exp[-2 KJz)(zb  -  z)]}/ 

{1  -  exp[— 2 Kd{z)(zb  -  z)]}.  (43) 

The  values  of  R  Y_  predicted  by  the  AM,  where  R  ,_  = 
R(z),  the  DMRM,  and  the  SMRM  are  shown  as  a  func¬ 
tion  of  depth  in  Figs.  14  and  15.  Estimates  of  R  Y  from 
the  SMRM  require  an  estimate  of  however,  it  was 
found  that  these  estimates  are  not  very  sensitive  to  Ky  , 
and  the  approximation  K,  =  KJj)  was  used  in  Figs.  14 
and  15.  It  can  be  seen  that  the  DMRM  method  does 
poorly  compared  with  the  AM  near  the  surface  but 
yields  much  better  estimates  than  the  AM  below  mid¬ 
depth.  The  SMRM  method  performs  better  than  the 
DMRM  near  the  surface,  but  performs  worse  than  the 
DMRM  near  the  bottom.  Both  the  DMRM  and  the 
SMRM  give  the  best  estimates  of  R,_  at  mid-depths, 
away  from  both  boundaries. 

Estimates  of  RA t)  from  the  DMRM  method  are 
more  susceptible  to  noise  in  the  irradiance  measure¬ 
ments  than  are  those  from  the  SMRM  or  AM.  How¬ 
ever,  smoothing  of  RA t)  was  found  to  be  effective  in 


R 


Fig.  16.  Estimation  of  R„  in  three  distinct  water  layers.  From 
top  to  bottom,  o)0  =  0.7,  0.65,  and  0.60.  Shown  are  R.Y  from  Eq. 
(13),  the  local  irradiance  ratio  R( t),  and  the  estimate  of-fiLA)  from 
the  DMRM. 

reducing  this  noise  in  simulated  data,  making  good 
estimates  of  R from  the  DMRM  possible. 

Profiles  of  KJz)  tend  to  be  less  influenced  by  the 
bottom  depth  and  albedo  than  are  R(z),  especially  if 
the  bottom  is  highly  absorbing.  Similar  to  R(z),  K(z) 
tends  to  be  closest  to  K,  at  its  minimum,  maximum, 
or  inflection  point.  Unfortunately,  a  shallow-water 
method  such  as  those  for  R  ,_  could  not  be  developed 
for  K,_. 

9.  IOP  Estimation  in  Inhomogeneous  Waters 

For  inhomogeneous  waters,  R,_  and  KY  as  computed 
from  Eqs.  (13)  and  (14)  now  depend  on  the  local  IOP’s 
at  z,  so  that  RJz)  and  KJz)  are  a  function  of  depth. 
If  the  optical  properties  vary  only  gradually  with 
depth,  then  R(z)  RJz)  and  KJz)  KJz)  below  the 
depths  where  the  surface  illumination  dominates. 
On  the  other  hand,  if  the  optical  properties  are  highly 
variable  with  depth,  some  type  of  weighted  average  of 
the  IOP’s  can  be  estimated,  but  the  fine  structure 
cannot  be  determined  accurately. 

Figure  16  shows  a  simulated  example  of  three  dis¬ 
tinct  water  layers.  For  all  three  layers  g  =  0.85,  but 
co0  =  0.7,  0.65,  and  0.60,  from  top  to  bottom.  Shown 
are  R  ,  from  Eq.  (13),  the  local  irradiance  ratio  R{ t), 
and  the  estimate  of  RJ t)  from  the  DMRM.  The  AM 
for  estimating  RAJ  in  the  three  layers  could  be  ap¬ 
plied  by  one  taking  the  maximum  R{ t)  in  the  upper 
layer,  the  inflection  point  in  the  middle  layer,  and  the 
asymptotic  value  in  the  deep  bottom  layer.  Because 
the  top  two  layers  in  this  case  are  both  affected  by  a 
distinctly  different  layer  below  them,  estimates  of  R  Y 
from  the  DMRM  are  more  accurate  than  those  from 
the  AM  just  above  the  interfaces,  but  are  less  accu¬ 
rate  just  below  the  interfaces.  Both  methods  do  a 
good  job  of  identifying  the  location  of  the  interfaces. 

10.  Summary 

We  have  tested  numerically  an  inverse  radiative 
transfer  method  for  estimating  two  inherent  optical 
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properties  from  the  irradiance  ratio  R(z)  and  the 
downward  diffuse  attenuation  coefficient  Kd(z)  com¬ 
puted  from  upward  and  downward  irradiance  mea¬ 
surements.  The  method  involves  two  steps.  First, 
the  IOP’s  R,  and  K,  are  estimated  from  R(z)  and 
Kd(z).  Second,  IOP’s  such  as  a  and  bb  are  calculated 
from  the  analytical  equations  in  (13)  and  (14)  for  A, 
and  K,_. 

The  simplest  approach  to  the  first  step  is  to  use  the 
AM,  where  it  is  assumed  R „  =  R(z)  and  K„  =  Kd(z). 
In  deep  waters,  the  exponential  method  of  Eqs.  (18) 
and  (19)  can  yield  better  results  than  the  AM  if  only 
shallow  irradiance  measurements  are  available.  On 
the  other  hand,  in  shallow  waters  the  DMRM  and 
SMRM  reflectance  models  of  Eqs.  (42)  and  (43)  yield 
better  estimates  for  R  ,  than  the  AM,  but  we  were 
unable  to  develop  a  deep-  or  shallow-measurement 
diffuse  attenuation  coefficient  model. 

In  the  second  step,  a  scattering  phase  function  is 
assumed,  and  w0  is  estimated  from  A,  with  an  iter¬ 
ative  search.  The  values  of  a  and  bb  are  then  calcu¬ 
lated  from  oj0  and  K,  .  This  approach  was  compared 
with  Haltrin’s  model  that  relates  a  and  bb  to  R  ,  and 
K,_  through  two  closed-form  equations  (39)  and  (41). 
Although  Haltrin’s  model  is  easier  to  implement,  it  is 
generally  less  accurate  because  no  a  priori  informa¬ 
tion  about  the  scattering  phase  function  can  be  incor¬ 
porated. 

Example  numerical  calculations  and  a  sensitivity 
analysis  have  shown  that  the  absorption  coefficient 
and  the  Gordon  parameter  G  =  bb/(a  +  bh)  can  be 
estimated  quite  well  even  if  the  scattering  asymme¬ 
try  factor  is  unknown.  However,  it  was  shown  with 
use  of  the  two-term  Henyey-Greenstein  phase  func¬ 
tion  that,  because  bb  is  sensitive  to  the  backscattering 
portion  of  the  phase  function,  it  is  important  to  use  a 
realistic  scattering  phase  model,  such  as  a  Petzold 
phase  function,  in  the  inverse  solution. 

Appendix  A:  Deep-Water  Asymptotic  Irradiance  Ratio 
and  Diffuse  Attenuation  Coefficient 

The  radiance  for  a  source-free  (i.e.,  no  fluorescence  or 
Raman  scattering  effects)  homogeneous  medium  can 
be  expressed  as  a  superposition  of  eigenmodes  of  Eq. 
(9).  The  downward  irradiance  Ed(z)  =  E+(z)  and  the 
upward  irradiance  Eu(z)  =  E_(z )  are  obtained  from  an 
integration  of  the  radiance  over  p,  with  Eqs.  (6)  and  (7): 

E±(z)  =  2  [C(vJ)gi(±vJ)exp(-cz/vj) 

j= i 

+  Ci-vjig^+vjiexpicz/vj)],  (Al) 


and  the  expansion  coefficients  C{±vj)  depend  on  the 
IOP’s  and  the  surface  illumination  and  bottom  albedo 
boundary  conditions.  Beyond  a  couple  of  optical 
depths  from  either  of  the  boundaries,  the  eigenvalues 
±vx  make  the  dominant  contributions  so  Eq.  (Al)  can 
be  approximated  by 

E±{z)  =  Ctvi^iCFviiexpi-cz/iq) 

+  C^-v/gZ+v/expicz /v/,  (A3) 

and  for  deep  waters  where  there  are  no  bottom  ef¬ 
fects,  C(-v j)  — >  0  and 

R,  =  gi(~v\)/gi(v/),  (A4) 

which  is  Eq.  (13).  In  a  similar  manner,  use  of  Eq. 
(A3)  with  Cf-iq)  =  0  and  Eq.  (8)  yields  of  Eq.  (14). 

Appendix  B:  Relationship  between  the  Sensitivity 
Coefficients  of  the  Forward  and  Inverse  Problems 

The  sensitivity  coefficients  for  the  inverse  problem 
(the  computation  of  co0  and  v1  from  R  ,  and  g )  can  be 
computed  more  easily  when  they  are  expressed  in 
terms  of  partial  derivatives  of  the  forward  problem 
(the  computation  of  R,_  and  v]  from  w0  and  g).  Let 
the  superscript  f  denote  forward-problem  variables, 
which  depend  on  oo0  and  g,  and  unmarked  variables 
denote  inverse-problem  variables.  If  we  make  a 
transformation  of  variables  and  use  the  Jacobian  of 
the  transformation,  it  follows  that 


5to0  1 

dA„  ~~  dRj/d O)0  ’ 


(Bl) 


dw0  _  1  dRj 

dg  ~  dRj/d co0  dg  ’ 


(B2) 


where  Rj is  computed  directly  from  Eq.  (13).  These 
results  are  consistent  with  those  of  Fig.  2.  Further¬ 
more,  from  the  chain  rule  and  the  previous  two  equa¬ 
tions, 


chi  dv/  du>0  dv/  1 

i)Ra=  du>0  ()R.r.  du>0  dRj/dM0  ’ 


(B3) 


Siq  dv/  dv/  dw0  dv/  dv/  1  dRj 
dg  dg  dw0  dg  dg  dto0  dRj / d«0  dg 


(B4) 


where 


g/vd  =  <Kv!,  |x)|xdp,, 

Jo 

^i(-vi)  =  |  cKvj,  p)|p|dp.  =  f  (K-Vj,  p.)pdp,  (A2) 

J-l  Jo 


Appendix  C:  Estimation  of  the  Asymptotic  Irradiance 
Ratio  Near  the  Bottom 

1 .  Derivation  of  the  DMRM 
From  Eqs.  (A3)  and  (A4), 

E+(z)  ±  EJz)^g1(v1)(  1  ±  RJ[C(v1)exp(-cz/v1) 

±  C(-V!)exp(c2/vi)].  (Cl) 


8696 


APPLIED  OPTICS  /  Vol.  36,  No.  33  /  20  November  1997 


After  squaring  the  last  equation  we  can  see  that 

[E+(z)  ±  E.{z)f  -^2(vi)(l  ±  R„)2[C2(Vl) 

X  exp(-2 cz/vj  4-  C^-iq) 

X  exp(2cz/vi)  ±  2C(vi)C(— vi)]. 

(C2) 


Equation  (42)  follows  after  we  evaluate  this  equation 
at  two  depths  z1  and  z2  and  take  the  difference  of  the 
results. 


2.  Derivation  of  the  SMRM 
From  Eqs.  (A3)  and  (A4), 


E_(z)  Ry,  +  C(— v1)exp(2c2v1)/C(v1) 
E+(z )  1  +  Ci—v^Ry,  exp(2cz/v1)/C(v1)  ’ 


(C3) 


so  it  follows  that 


R(Z)  ~  Ryy  = 


C(—  V!))!  —  Ry!)exp(2cz/v1)/C(v1) 


1  +  RXC(— v1)exp(2c2/r1)/C(r1) 
For  z  in  Eq.  (C4)  just  below  the  surface, 


(C4) 


R(  0  +  )  -  Ryy 


C(-Vl){l  -  RJ/CivJ 
1  +  RjC(~V  l)/C(Vl)  ’ 


(C5) 


whereas  for  2  in  Eq.  (C4)  at  the  bottom,  z  =  zb,  the 
reflectance  equals  the  bottom  albedo  Rb,  and 


_  C(-Tq)(l  -  i?Jexp(2cz6/v1)/C(v1) 
6  ”  1  +  RyyC(— v1)exp(2czi/r1)/C(v1) 

From  Eqs.  (C5)  and  (C6), 


(Rb  ~  R^)exp(-2czb/v1)  = 

[Jg(0+)-ig,][l  +  C(-v1)ig,/C(v1)] 

1  +.R„C(— v1)exp(2cz6/v1)/C(v1) 

IfRJ^-vJ/COq)  «  1,  then 

(Rb  -  Ryy)exp(-2czb/v1)  ~  R( 0+)  -  (C8) 

It  follows  that  if  the  bottom  albedo  is  known  and 
=  c/vr  is  estimated  from  Ed[z),  then  R„  can  be  esti¬ 
mated  from  R(  0+): 


R-yy  =  [R{ 0+)  -  Rb  exp(— 2if„z6)]/[l  -  exp(-2RLz6)]. 

(C9) 

This  equation  has  been  used  by  others,  e.g.,  Ref.  27. 
Equation  (43)  is  a  generalization  away  from  z  =  0+ 
and  uses  either  the  local  Kd(z)  or  the  best  estimate  of 

Kyy. 

This  research  was  supported  by  the  U.S.  Office  of 
Naval  Research  and  by  the  San  Diego  Supercomputer 
Center.  We  thank  Howard  Gordon,  Vladimir  Hal- 
trin,  Emmanuel  Boss,  Lydia  Sundman,  and  the  re¬ 
viewers  for  helpful  suggestions  and  Knut  Stamnes  for 
providing  us  with  the  computer  code  disortb. 


References 

1.  R.  W.  Preisendorfer,  Hydrologic  Optics,  V.l.  NTIS  PB  259793/ 
8ST  (National  Technical  Information  Service,  Springfield,  Va., 
1976). 

2.  D.  A.  Kiefer  and  B.  G.  Mitchell,  “A  simple,  steady-state  de¬ 
scription  of  phytoplankton  growth  based  on  absorption  cross 
section  and  quantum  efficiency,”  Limnol.  Oceanogr.  28,  770- 
776  (1983). 

3.  M.  Kishino,  “Interrelationship  between  light  and  phytoplank¬ 
ton  in  the  sea,”  in  Ocean  Optics,  R.  W.  Spinrad,  K.  L.  Carder, 
and  M.  J.  Perry,  eds.  (Oxford  U.  Press,  New  York,  1994). 

4.  A.  A.  Gershun,  “The  light  field,”  J.  Math.  Phys.  (Cambridge, 
Mass.)  18,  51-151  (1939). 

5.  K.  L.  Carder,  D.  J.  Collins,  M.  J.  Perry,  H.  L.  Clark,  J.  M. 
Mesias,  J.  S.  Cleveland,  and  J.  Greenier,  “The  interaction  of 
light  with  phytoplankton  in  the  marine  environment,”  in 
Ocean  Optics  VIII,  M.  A.  Blizard,  ed.,  Proc.  SPIE  637,  42-55 
(1986). 

6.  L.  Prieur  and  S.  Sathyendranath,  “An  optical  classification  of 
coastal  and  oceanic  waters  based  on  the  specific  spectral  ab¬ 
sorption  curves  of  phytoplankton  pigments,  dissolved  organic 
matter,  and  other  particulate  materials,”  Limnol.  Oceanogr. 
26,  671-689  (1981). 

7.  J.  R.  V.  Zaneveld,  “A  reflecting  tube  absorption  meter,”  in 
Ocean  Optics X,  R.  W.  Spinrad,  ed.,  Proc.  SPIE  1302,  124-136 
(1990). 

8.  H.  R.  Gordon,  “Modeling  and  simulating  radiative  transfer  in 
the  ocean,”  in  Ocean  Optics,  R.  W.  Spinrad,  K.  L.  Carder,  and 
M.  J.  Perry,  eds.  (Oxford  U.  Press,  New  York,  1994). 

9.  H.  R.  Gordon  and  G.  C.  Boynton,  “A  radiance-irradiance  in¬ 
version  algorithm  for  estimating  the  absorption  and  backscat- 
tering  coefficients  of  natural  waters:  homogeneous  waters,” 
Appl.  Opt.  36,  2636-2641  (1997). 

10.  Z.  Tao,  N.  J.  McCormick,  and  R.  Sanchez,  “Ocean  source  and 
optical  property  estimation  using  explicit  and  implicit  algo¬ 
rithms,”  Appl.  Opt.  33,  3265-3275  (1994). 

11.  S.  Chandrasekhar,  Radiative  Transfer  (Oxford  U.  Press,  New 
York,  1950). 

12.  M.  Benassi,  R.  D.  M.  Garcia,  A.  H.  J.  Karp,  and  C.  E.  Siewert, 
“A  high-order  spherical  harmonics  solution  to  the  standard 
problem  in  radiative  transfer,”  Astrophys.  J.  280,  853-864 
(1984). 

13.  N.  J.  McCormick,  “Asymptotic  optical  attenuation,”  Limnol. 
Oceanogr.  37,  1570-1578  (1992). 

14.  C.  E.  Siewert,  “The  FN  method  for  solving  radiative-transfer 
problems  in  plane  geometry,”  Astrophys.  Space  Sci.  58,  131- 
137  (1978). 

15.  L.  C.  Henyey  and  J.  L.  Greenstein,  “Diffuse  radiation  in  the 
galaxy,”  Astrophys.  J.  93,  70-83  (1941). 

16.  K.  Stamnes,  “The  Chandrasekhar  method  and  its  applications 
to  atmospheric  radiative  transfer,”  Trans.  Am.  Nucl.  Soc.  71, 
213-214  (1994). 

17.  Z.  Jin  and  K.  Stamnes,  “Radiative  transfer  in  nonuniformly 
refracting  layered  media  such  as  the  atmosphere/ocean  sys¬ 
tem,”  Appl.  Opt.  33,  431-442  (1994). 

18.  N.  J.  McCormick,  “Analytical  transport  theory  applications  in 
optical  oceanography,”  Ann.  Nucl.  Energy  23,  381-395  (1996). 

19.  J.  R.  V.  Zaneveld,  “An  asymptotic  closure  theory  for  irradiance 
in  the  sea  and  its  inversion  to  obtain  the  inherent  optical 
properties,”  Limnol.  Oceanogr.  34,  1442-1452  (1989). 

20.  N.  J.  McCormick,  “Mathematical  models  for  the  mean  cosine  of 
irradiance  and  the  diffuse  attenuation  coefficient,”  Limnol. 
Oceanogr.  40,  1013-1018  (1995). 

21.  H.  R.  Gordon,  O.  B.  Brown,  and  M.  M.  Jacobs,  “Computed 
relationships  between  the  inherent  and  apparent  optical  prop¬ 
erties  of  a  flat  homogeneous  ocean,”  Appl.  Opt.  14,  417-427 
(1975). 

22.  W.  H.  Press,  B.  P.  Flannery,  S.  A.  Teukolsky,  and  W.  T.  Vet- 


20  November  1997  /  Vol.  36,  No.  33  /  APPLIED  OPTICS 


8697 


terling,  Numerical  Recipes  (Cambridge  U.  Press,  New  York, 
1989),  pp.  274-286. 

23.  T.  J.  Petzold,  “Volume  scattering  functions  for  selected  ocean 
waters,”  SIO  Ref.  72-78  (Scripps  Institution  of  Oceanography, 
Visibility  Laboratory,  San  Diego,  Calif.,  1972). 

24.  P.  W.  Francisco  and  N.  J.  McCormick,  “Chlorophyll  concentra¬ 
tion  effects  on  asymptotic  optical  attenuation,”  Limnol.  Ocean- 
ogr.  39,  1195-1205  (1994). 

25.  G.  W.  Kattawar,  “A  three-parameter  analytic  phase  function 
for  multiple  scattering  calculations,”  J.  Quant.  Spectrosc.  Ra- 
diat.  Transfer  15,  839-849  (1975). 


26.  V.  I.  Haltrin,  “Algorithm  for  computing  apparent  optical  prop¬ 
erties  of  shallow  waters  under  arbitrary  surface  illumination,” 
in  Proceedings  of  the  International  Airborne  Remote  Sensing 
Conference  and  Exhibition  (Environmental  Research  Institute 
of  Michigan,  Ann  Arbor,  Mich.,  1997),  pp.  463-470. 

27.  N.  T.  O’Neill,  A.  R.  Kalinauskas,  J.  D.  Dunlop,  A.  B.  Hollinger, 
H.  Edel,  M.  Casey,  and  J.  Gibson,  “Bathymetric  analysis  of 
geometrically  corrected  imagery  data  collected  using  a  two 
dimensional  imager,”  in  Ocean  Optics  VIII,  M.  A.  Blizard,  ed., 
Proc.  SPIE  637,  191-202  (1986). 


8698 


APPLIED  OPTICS  /  Vol.  36,  No.  33  /  20  November  1997 


